function verify
close all
global delta Nx

delta = 1e-2;

L = @(tau) diag([1 - 2 * tau; -2 * tau; -1 - 2 * tau]);
% Q = @(tau) [1 -1 0; 2 1 -1; 0 -1 -1];

A = @(tau) [tau - 1 -1 0; 2 tau -1; 0 -1 -1];
A = @(tau) inv(Q(tau)) * L(tau) * Q(tau);
U = [-1; 1; 1];


A = @(tau)[1-tau,1;1,-1-tau];
U = [1; 1];

Nx = size(A(0), 1);
Nmodes = Nx;


% tau = [0 logspace(-5, 0, 200)];
tau = linspace(0, 1, 500);
N = length(tau);

dt = min(diff(tau));

% the initial condition


opts = odeset('RelTol',1e-12);
num = ode15s(@(t, x) A(delta * t) * x, [0, tau(end) / delta], U, opts);


% compute eigenvalues, eigenfunctions, adjoints
lambda = zeros(Nmodes, N);
lambda_a = lambda;
v = zeros(Nx, Nmodes, N);
adj = v;
for i = 1:N
    
    
    
    [ef, ev] = eigs(A(tau(i)));
    [tmp, P] = sort(diag(ev), 'descend');
    lambda(:,i) = tmp;
    v(:,:,i) = ef(:,P);
    
    [ef, ev] = eigs(A(tau(i))');
    [tmp, P] = sort(diag(ev), 'descend');
    lambda_a(:,i) = tmp;
    adj(:,:,i) = ef(:,P);

    
    % normalize eigenfunctions
    for j = 1:Nmodes
        v(:,j,i) = v(:,j,i) / sqrt(ip(v(:,j,i), v(:,j,i)));
        adj(:,j,i) = adj(:,j,i) / ip(adj(:,j,i), v(:,j,i));
    end
    
end

% go through and correct sign flips in the eigenfunctions
for i = 2:N
    for j = 1:Nmodes
        if (v(1,j,i) * v(1,j,i-1) < 0 && v(Nx,j,i) * v(Nx,j,i-1) < 0)
            v(:,j,i) = -v(:,j,i);
        end
        
        if (adj(1,j,i) * adj(1,j,i-1) < 0 && adj(Nx,j,i) * adj(Nx,j,i-1) < 0)
            adj(:,j,i) = -adj(:,j,i);
        end
        
    end
end
    
    

% interpolate eigenvectors for easy differentiation
v_i = @(t_i) interp3(1:Nmodes, 1:Nx, tau, v, (1:Nmodes)', (1:Nx), t_i);


% plot the eigenvalues
figure
subplot(2,1,1);
plot(tau, lambda);
xlabel('\tau');
ylabel('Eigenvalues');
subplot(2,1,2);
plot(tau, lambda_a);
xlabel('tau');
ylabel('Adjoint eigenvalues');

% extract coeffs
num_i = deval(num, tau / delta);
c_num = zeros(Nmodes, N);
for i = 1:N
    for j = 1:Nmodes
        c_num(j,i) = ip(reshape(adj(:,j,i), Nx, 1), num_i(:,i));
    end
end

% compute the c_i
c_i = zeros(Nmodes, 1);
for i = 1:Nmodes
    c_i(i) = ip(reshape(adj(:,i,1), Nx, 1), U);
end
fprintf('c_i:\n');
disp(c_i);


% loop again to compute gamma
gamma = zeros(Nmodes, Nmodes, N);
for i = 1:N
   
    if (i == 1)
        dv = (v_i(tau(i) + dt) - v_i(tau(i))) / dt;
    elseif (i == N)
        dv = (v_i(tau(i)) - v_i(tau(i) - dt)) / dt;
    else
        dv = (v_i(tau(i) + dt) - v_i(tau(i) - dt)) / 2 / dt;
    end
    
    % gamma(j,k,.) = v*(j)' * v(k)
    
    for j = 1:Nmodes
        for k = 1:Nmodes
            gamma(j,k,i) = -ip(reshape(adj(:,j,i), Nx, 1), dv(:,k));
        end
    end
    
end

% plot gamma
figure
for j = 1:Nx
    for k = 1:Nx
        subplot(Nx, Nx, k + (j - 1) * Nx);
        plot(tau, reshape(gamma(j,k,:), N, 1));
        xlabel('\tau');
        ylabel(sprintf('\\gamma_{%d%d}',j,k));
    end
end

% compute beta
tmp = zeros(N, Nmodes);
for i = 1:Nmodes
    tmp(:,i) = reshape(gamma(i,i,:), N, 1);
end
beta = cumtrapz(tau', lambda' / delta + tmp)';

% plot beta
figure
plot(tau, beta);
xlabel('\tau');
ylabel('\beta');

% compute leading order coeffs
c_lo = zeros(Nmodes, N);
for i = 1:Nmodes 
    c_lo(i,:) = c_i(i) * exp(beta(i,:));
end

fprintf('c_lo(t = 0):\n');
disp(c_lo(:,1));
fprintf('c_num(t = 0):\n');
disp(c_num(:,1));

% do the corrections
c_corr = zeros(Nmodes, N);
for j = 1:Nmodes
    for l = 1:Nmodes
        if (l ~= j)
            int = cumtrapz(tau, reshape(gamma(j, l,:) .* gamma(l, j,:), 1, N) ./ (lambda(j,:) - lambda(l,:)));
            
            tmp = -c_i(l) * reshape(gamma(j, l, :), 1, N) ./ (lambda(j, :) - lambda(l, :)) .* exp(beta(l,:)) + ...
                (c_i(j) * int + gamma(j, l, 1) * c_i(l)  / (lambda(j,1) - lambda(l,1))) .* exp(beta(j,:));
                        
            c_corr(j,:) = c_corr(j,:) + tmp;
        end
    end
end

% plot the solutions
A_num = num_i;
A_lo = zeros(Nx,N);
A_corr = A_lo;
for i = 1:N
    A_lo(:,i) = reshape(v(:,:,i), Nx, Nmodes) * c_lo(:,i);
    A_corr(:,i) = reshape(v(:,:,i), Nx, Nmodes) * (c_lo(:,i) + delta * c_corr(:,i));
end
figure
for j = 1:Nmodes
    subplot(ceil(Nmodes/2),2,j);
    plot(tau, A_num(j,:), 'k', tau, A_lo(j,:), 'b', tau, A_lo(j,:) + delta * A_corr(j,:), 'r');
    xlabel('\tau');
    ylabel(strcat('u_',num2str(j)));
end

% plot the coefficients
figure
for j = 1:Nmodes
    subplot(ceil(Nmodes/2),2,j);
    plot(tau, c_num(j,:), 'k', tau, c_lo(j,:), 'b', tau, c_lo(j,:) + delta * c_corr(j,:), 'r');
    xlabel('\tau');
    ylabel(strcat('c_',num2str(j)));
end

% plot the error
figure
for j = 1:Nmodes
    subplot(ceil(Nmodes/2),2,j);
    tmp = c_lo(j,:) + delta * c_corr(j,:);
    err = abs((tmp - c_num(j,:)) ./ c_num(j,:));
    semilogy(tau, err);
    xlabel('\tau');
    ylabel(strcat('err_',num2str(j)));
end


%--------------------------------------------------------------------------

function v = ip(f1, f2)

v = f1' * f2;